********************************************************************************
*** Data Paper (MIPD Release 2.0): State analysis
***
*** Created: 3-17-23
*** Modified:
***
********************************************************************************

*** Make sure to install spmap and shp2dta
*ssc install spmap
*ssc install shp2dta

*** Load the data
use "State Disproportionality.dta", clear

********************************************************************************
*** Display the biggest categories across states
********************************************************************************

* Display the most common specific categories in each state, ranked.
preserve
	drop if inrange(mipcode, 1301, 1599)

	duplicates drop state mipcode, force
	bys state: egen rank = rank(perc), field

	sort rank
	levelsof state, local(STATE)
	foreach s of local STATE {
		list state mipcode perc rank if state == `s' & inrange(rank, 1, 5)
	}
	
	bys state: egen rank_disp = rank(disp), field

	sort rank_disp
	levelsof state, local(STATE)
	foreach s of local STATE {
		list state mipcode perc p disp rank_disp if state == `s' & inrange(rank_disp, 1, 5)
	}	
restore


********************************************************************************
*** Spatial dependence in problem disproportionality across states
********************************************************************************

* Simple 48x48 matrix (excludes Alaska, DC and Hawaii)
spatwmat using "Maps\Wstates.dta", name(W)

postfile spat mipcode geary geary_p moran moran_p using "Spatial Clustering.dta", replace

*** Spatial association across specific codes
levelsof mipcode if !inlist(mipcode, 206, 213, 214, 215, 216, 301, 302, 306, 307, 402, 406, 609, 701, 704, 705, 802, 904, 906, 908, 909, 1009, 1013, 1104, 1105, 1106, 1201, 1301, 1302, 1401, 1501, 1502, 1503, 1504), local(MIP)
qui foreach mip of local MIP {
	preserve
		nois display "MIP code = `mip'" 
		
		keep if mipcode == `mip'
	
		drop if inlist(state, 2, 11, 15)
		duplicates drop state, force
		sort state
	
		spatgsa perc, w(W) m g two
		mat G = r(Geary)
		mat M = r(Moran)
		
		local geary = G[1,1]
		local geary_p = G[1,5]
		local moran = M[1,1]
		local moran_p = M[1,5]
		
		post spat (`mip') (`geary') (`geary_p') (`moran') (`moran_p')
	restore
}

postclose spat

preserve
	use "Spatial Clustering.dta", clear
	
	sort moran
	list mipcode geary geary_p moran moran_p
	
	* What percentage have *positive* Moran's I?
	gen positive = cond(moran > 0, 1, 0)
	tab positive
	
	* What percentage have *statistically significant* Moran's I?
	gen ss = cond(moran_p <= 0.05, 1, 0)
	tab ss
	
	* Moran's I for Racism
	list moran moran_p if mipcode == 309
	
	* Moran's I for Immigration
	list moran moran_p if mipcode == 509
restore

********************************************************************************
*** Spatial map of problem disproportionality across states
********************************************************************************
shp2dta using "Maps\cb_2018_us_state_500k", database(us) coordinates(uscoord) genid(id) replace 

*** Open up the us.dta file and describe
preserve
	use us.dta, clear
	desc

	destring STATEFP, force gen(state)
	
	drop if inlist(state, 2, 15, 60, 66, 69, 72, 78)
	sort state

	tempfile us
	save `us', replace 
restore	
	
*** Open up the state disproportionaliy file
use "State Disproportionality.dta", clear

duplicates drop state mipcode, force

sort state
merge state using `us'
drop if _merge == 2
drop _merge

*** Map out problem disproportionality of various issue categories across states
spmap disp using uscoord if mipcode == 309, id(id) fcolor(Reds) title("Racism/Discrimination")
spmap disp using uscoord if mipcode == 509, id(id) fcolor(Reds) title("Immigration")

*** Maps in the Online Appendix
spmap disp using uscoord if mipcode == 903, id(id) fcolor(Reds) title("Religion")
spmap disp using uscoord if mipcode == 508, id(id) fcolor(Reds) title("Agricultural")
spmap disp using uscoord if mipcode == 102, id(id) fcolor(Reds) title("Employment")
spmap disp using uscoord if mipcode == 801, id(id) fcolor(Reds) title("Environment")
